Forest Management Study¶
Clark, C. M., Richkus, J., Jones, P. W., Phelan, J., Burns, D. A., de Vries, W., Du, E., Fenn, M. E., Jones, L., & Watmough, S. A. (2019). A synthesis of ecosystem management strategies for forests in the face of chronic nitrogen deposition. Environmental Pollution, 248, 1046–1058. https://doi-org.proxy.lib.miamioh.edu/10.1016/j.envpol.2019.02.006¶
In [1]:
# You need to run this block of code to import the libraries that are needed for the data manipulations, visualizations, etc.
from pydoc import help # can type in the python console `help(name of function)` to get the documentation
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import scale
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from scipy import stats
from IPython.display import display, HTML
Import data, data prep, and EDA¶
In [2]:
# Read csv file into pandas dataframe
forest = pd.read_csv("ForestData.csv")
In [3]:
# Examine first few cases
forest.head()
Out[3]:
In [4]:
# List the columns
list(forest.columns)
Out[4]:
In [5]:
forest.dtypes
Out[5]:
Filtering out non relevant cases¶
In [6]:
df = forest[forest['Ecosystem.element.assessed..soil.plant.microbe.fungi.fauna.'].isin(['Plant']) ]
In [7]:
df
Out[7]:
In [8]:
a=['non-native' , 'alien species']
df = df[~df['Endpoint.Modifier'].isin(a)]
In [9]:
df
Out[9]:
In [10]:
df = df[df['Endpoint'].isin(['shannon diversity', 'simpson diversity', 'species richness']) ]
In [11]:
df
Out[11]:
Removing missing values¶
In [12]:
# Getting rid of rows with missing data on the response variable
df = df[df['Response..Measurement'].notna()]
In [13]:
# Getting rid of rows with missing data on the response variable standard deviation
df = df[df['Response.Standard.Deviation'].notna()]
In [14]:
# Getting rid of rows with missing data on the response variable sample size
df = df[df['Response.Number.of.Observations'].notna()]
In [15]:
# Getting rid of rows with missing data on the control variable standard deviation
df = df[df['Control.Reference.Measurement'].notna()]
In [16]:
# Getting rid of rows with missing data on the response variable standard deviation
df = df[df['Control.Standard.Deviation'].notna()]
In [17]:
# Getting rid of rows with missing data on the response variable standard deviation
df = df[df['Control.Number.of.Observations'].notna()]
In [18]:
df
Out[18]:
Parsing out year from reference to remove prior to 1996¶
In [19]:
df['Year'] = df['Reference'].astype('str').str.extractall('(\d+)').unstack().fillna('').sum(axis=1).astype(int)
In [20]:
df['Year'] = df['Year'].astype(float).apply(lambda x:x if (x > 0) else 9999)
In [21]:
df[df['Year'] > 1995]
Out[21]:
In [22]:
# Create a function that we can re-use
# This block of code will create a histogram. You don't have to use all of this code next time you want to create one
def show_distribution(var_data):
from matplotlib import pyplot as plt
# Get statistics
min_val = var_data.min()
max_val = var_data.max()
mean_val = var_data.mean()
med_val = var_data.median()
mod_val = var_data.mode()[0]
print('Minimum:{:.2f}\nMean:{:.2f}\nMedian:{:.2f}\nMode:{:.2f}\nMaximum:{:.2f}\n'.format(min_val,
mean_val,
med_val,
mod_val,
max_val))
# Create a figure for 2 subplots (2 rows, 1 column)
fig, ax = plt.subplots(2, 1, figsize = (10,4))
# Plot the histogram
ax[0].hist(var_data)
ax[0].set_ylabel('Frequency')
# Add lines for the mean, median, and mode
ax[0].axvline(x=min_val, color = 'gray', linestyle='dashed', linewidth = 2)
ax[0].axvline(x=mean_val, color = 'cyan', linestyle='dashed', linewidth = 2)
ax[0].axvline(x=med_val, color = 'red', linestyle='dashed', linewidth = 2)
ax[0].axvline(x=mod_val, color = 'yellow', linestyle='dashed', linewidth = 2)
ax[0].axvline(x=max_val, color = 'gray', linestyle='dashed', linewidth = 2)
# Plot the boxplot
ax[1].boxplot(var_data, vert=False)
ax[1].set_xlabel('Value')
# Add a title to the Figure
fig.suptitle('Data Distribution')
# Show the figure
fig.show()
col = df['Response..Measurement']
# Call the function
show_distribution(col)
In [23]:
# Get the variable to examine
col = df['Control.Reference.Measurement']
# Call the function
show_distribution(col)
In [24]:
# Exploring relationships - this code will create scatter plots
import matplotlib.pyplot as plt
from pandas.plotting import scatter_matrix
# selecting numerical features
features = ['Response..Measurement', 'Control.Reference.Measurement', 'Response.Standard.Deviation', 'Control.Standard.Deviation']
# plotting the scatter matrix
scatter_matrix(df[features], figsize=(12,12))
plt.xticks(rotation=90)
plt.show()
In [25]:
# If you want to compare groups based on one or more quantitative variables, use this code
col_list = ['Management.Method', 'Response..Measurement', 'Control.Reference.Measurement']
short_df = df[col_list]
import seaborn as sns
rs=1999
df_long = pd.melt(short_df.sample(150,random_state=rs), "Management.Method", var_name="Columns", value_name="Values")
f,ax = plt.subplots(figsize=(16,8))
#plt.xticks(rotation=90)
plt.ylim(0, 125)
#plt.xlim(0, None)
sns.boxplot(x="Columns", y="Values", hue="Management.Method", data=df_long)
Out[25]:
In [26]:
df['ts2'] = df['Response.Standard.Deviation']**2
In [27]:
df['cs2'] = df['Control.Standard.Deviation']**2
In [28]:
df['term1'] = ((df['Response.Number.of.Observations'] - 1) * (df['ts2']))
In [29]:
df['term2'] = ((df['Control.Number.of.Observations'] - 1) * (df['cs2']))
In [30]:
df['numerator'] = df['term1'] + df['term2']
In [31]:
df['denominator'] = df['Response.Number.of.Observations'] + df['Control.Number.of.Observations'] - 2
In [32]:
df['divided'] = df['numerator'] / df['denominator']
In [33]:
df['pooledVar']=df['divided']**(1/2)
In [34]:
df['target'] = ((df['Response..Measurement'] - df['Control.Reference.Measurement']) / (df['pooledVar']))
In [35]:
# Getting rid of rows with missing data on the response variable
df = df[df['target'].notna()]
In [36]:
df
Out[36]:
In [37]:
# If you want to compare groups based on one or more quantitative variables, use this code
col_list = ['Management.Method', 'target']
short_df = df[col_list]
import seaborn as sns
rs=1999
df_long = pd.melt(short_df.sample(150,random_state=rs), "Management.Method", var_name="Management Method", value_name="Effect Size")
f,ax = plt.subplots(figsize=(16,8))
#plt.xticks(rotation=90)
plt.ylim(-10, 10)
#plt.xlim(0, None)
sns.boxplot(x="Management Method", y="Effect Size", hue="Management.Method", data=df_long)
Out[37]:
In [38]:
# If you want to compare groups based on one or more quantitative variables, use this code
col_list = ['Ecosystem.element.assessed..soil.plant.microbe.fungi.fauna.', 'target']
short_df = df[col_list]
import seaborn as sns
rs=1999
df_long = pd.melt(short_df.sample(150,random_state=rs), "Ecosystem.element.assessed..soil.plant.microbe.fungi.fauna.", var_name="Columns", value_name="Values")
f,ax = plt.subplots(figsize=(16,8))
#plt.xticks(rotation=90)
plt.ylim(-10, 10)
#plt.xlim(0, None)
sns.boxplot(x="Columns", y="Values", hue="Ecosystem.element.assessed..soil.plant.microbe.fungi.fauna.", data=df_long)
Out[38]:
In [39]:
#Keeping only the columns we want
cols_to_keep = ['Reference',
'Habitat.Type',
'Management.Method',
'Management.Effect.Type',
'Ecosystem.element.assessed..soil.plant.microbe.fungi.fauna.',
'Endpoint',
'Measurement.Unit',
'Endpoint.Modifier',
'Magnitude.or.degree.of.change',
'Magnitude.or.degree.of.change..text.',
'Response..Measurement',
'Response.Number.of.Observations',
'Response.Standard.Deviation',
'Control.Reference.Measurement',
'Control.Number.of.Observations',
'Control.Standard.Deviation',
'general.area',
'Location',
'Year',
'target']
df2=df[cols_to_keep]
In [40]:
#Create dummy variables
combine = pd.concat([df2,pd.get_dummies(df['Management.Method'], prefix='ID')],axis=1)
combine.columns
Out[40]:
In [41]:
#Create dummy variables
combine2 = pd.concat([combine,pd.get_dummies(combine['Habitat.Type'], prefix='ID')],axis=1)
combine2.columns
Out[41]:
In [42]:
combine3 = pd.concat([combine2,pd.get_dummies(combine2['Management.Effect.Type'], prefix='ID')],axis=1)
combine3.columns
Out[42]:
In [43]:
combine4 = pd.concat([combine3,pd.get_dummies(combine3['Location'], prefix='ID')],axis=1)
combine4.columns
Out[43]:
In [44]:
combine5 = pd.concat([combine4,pd.get_dummies(combine4['Endpoint'], prefix='ID')],axis=1)
combine5.columns
Out[44]:
Decision trees¶
In [45]:
combine5.head(10)
Out[45]:
In [46]:
#conda install graphviz
In [47]:
#conda update -n base -c defaults conda
In [48]:
# Creating a list of features for the decision tree
# Include in this list, all the predictor variables that you want to use in your decision tree
features=['ID_Carbon Addition',
'ID_Prescribed Burn',
'ID_Thinning',
'ID_Coniferous Forest',
'ID_Deciduous Forest',
'ID_Forest',
'ID_community composition',
'ID_shannon diversity',
'ID_simpson diversity',
'ID_species richness',
]
In [49]:
target=['target']
In [50]:
# This block of code runs the decision tree analysis
# Replace TargetVar with the name of the target or dependent variable of interest
# max_depth and max_leaf_nodes are hyper parameters that are set by the researcher, they limit the size of the tree
from sklearn.tree import DecisionTreeRegressor, export_graphviz
from sklearn import datasets, tree
y = combine5["target"]
X = combine5[features]#.astype(float)
dt = DecisionTreeRegressor(min_samples_split=10, min_samples_leaf=10, random_state=99, max_depth=5, max_leaf_nodes=7)
dt.fit(X, y)
Out[50]:
In [51]:
text_representation = tree.export_text(dt)
print(text_representation)
In [52]:
with open("decistion_tree.log", "w") as fout:
fout.write(text_representation)
In [53]:
fig = plt.figure(figsize=(25,20))
_ = tree.plot_tree(dt,
feature_names=features,
class_names=target,
filled=True)
In [54]:
# If you want to compare groups based on one or more quantitative variables, use this code
col_list = ['Habitat.Type', 'target']
short_df = df[col_list]
import seaborn as sns
rs=1999
df_long = pd.melt(short_df.sample(150,random_state=rs), "Habitat.Type", var_name="Habitat Type", value_name="Effect Size")
f,ax = plt.subplots(figsize=(16,8))
#plt.xticks(rotation=90)
plt.ylim(-10, 10)
#plt.xlim(0, None)
sns.boxplot(x="Habitat Type", y="Effect Size", hue="Habitat.Type", data=df_long)
Out[54]:
In [55]:
#conda install -c plotly plotly=5.10.0
In [56]:
# selecting rows based on condition
richness = df[df['Endpoint'] == 'species richness']
In [57]:
richness.rename(columns = {'target': 'Effect Size', 'Management.Method':'Management Method', 'Habitat.Type':'Habitat Type'}, inplace = True)
richness
Out[57]:
In [61]:
import plotly.express as px
In [62]:
px.box(richness, x="Effect Size", y="Management Method", orientation="h", color="Habitat Type", notched=True,
category_orders={"Management Method": ["Carbon Addition", "Prescribed Burn", "Thinning"]})
In [63]:
# selecting rows based on condition
shannon = df[df['Endpoint'] == 'shannon diversity']
In [64]:
shannon.rename(columns = {'target': 'Effect Size', 'Management.Method':'Management Method', 'Habitat.Type':'Habitat Type'}, inplace = True)
richness
Out[64]:
In [65]:
px.box(shannon, x="Effect Size", y="Management Method", orientation="h", color="Habitat Type", notched=True,
category_orders={"Management Method": ["Carbon Addition", "Prescribed Burn", "Thinning"]})
In [66]:
# selecting rows based on condition
simpson = df[df['Endpoint'] == 'simpson diversity']
In [67]:
simpson.rename(columns = {'target': 'Effect Size', 'Management.Method':'Management Method', 'Habitat.Type':'Habitat Type'}, inplace = True)
richness
Out[67]:
In [68]:
px.box(simpson, x="Effect Size", y="Management Method", orientation="h", color="Habitat Type", notched=True,
category_orders={"Management Method": ["Carbon Addition", "Prescribed Burn", "Thinning"]})